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Abstract 


A transient finite-difference wave envelope 
formulation is presented for sound propagation, 
without steady flow. Before the finite difference 
equations are formulated, the governing wave equa- 
tion is first transformed to a form whose solution 
tends not to oscillate along the propagation direc- 
tion. This transformation reduces the required 
number of grid points by an order of magnitude. 
Physically, the transformed pressure represents 
the amplitude of the conventional sound wave. The 
derivation for the wave envelope transient wave 
equation and appropriate boundary conditions are 
presented as well as the difference equations and 
stability requirements. To illustrate the method, 
example solutions are presented for sound propaga- 
tion in a straight hard-wall duct and in a two- 
dimensional straight soft-wall duct. The numerical 
results are in good agreement with exact analytical 
results. 


Nomenclature 


a m ,b m’ 

d m’ e m’ 


c m , cell coefficients (table I] 
f m> 


m’ ''m’^m* 


9m»J 


l m ,rn m 
m m 

c o* 

adB 

E x 

f* 

H* 

I 

! x 

i 

J 

L* 

P 

P 

T* 

t 

at 

U 

u 

V 

x 

ax 

y 

Ay 


ambient speed of sound, m/sec 

decrease in decibels 

acoustic power 

frequency, Hz 

duct height, m 

number of axial grid points 

time averaged axial intensity 

VT 

number of transverse grid points 
length of duct, m . 
dimensionless pressure ^P*/p Q C*^ 
dimensionless wave envelope pressure, 
P(x,y,t) 
period, 1/f*, s 
dimensionless time, t*/T* 
time step 

dimensionless axial particle velocity 
wave envelope axial particle velocity 
dimensionless transverse particle 
velocity 

dimensionless axial coordinate, x*/H* 
axial grid spacing 

dimensionless traverse coordinate, y*/H* 
transverse grid spacing 


o 

Z* impedance, kg/m s 

0 coefficient (Eq. (25)) 

B coefficient (Eq. (26)) 

1 coefficient (Eq. (27)) 

c specific acoustic impedance 

n dimensionless frequency (Eq. (2)) 

n wave envelope frequency 

x effective wave length 

e dimensionless resistance 

* 3 

p Q ambient. air density, kg/m . 

x dimensionless reactance (for e +i “* t *) 

angular frequency, rad/S 

Subscripts 

e exit condition 

i axial grid index 

j transverse grid index 

m cell index 

max maximum stable time increment (Eq. (29)) 

o ambient condition 

Superscripts 

* dimensional quantity 

k time step 


Introduction 


"Steady-state" finite difference and finite 
element theories*^ have been developed to study 
sound propagation in free space and in complex 
ducts with axial variations in cross-sectional 
area, wall liner impedance (absorbers) and with 
gradients in flow Mach number. In the "steady 
state" theory, the pressure and acoustic veloci- 
ties are assumed to be simple harmonic functions 
of time; thus, the equations governing sound propa- 
gation (linearized gas dynamic equations’ (pg. 5) 
become independent of time. Generally the "steady 
state" finite difference and finite element numeri- 
cal algorithms have been limited to low frequencies 
and short ducts, because many grid points or ele- 
ments were required to resolve the axial wave 
length of sound and because of the large matrices 
associated with the solution of time independent 
equations. 

In many practical situations something (a sup- 
pressor’s impedance or a turboprop’s blade geome- 
try) needs to be optimized in some manner to obtain 
the maximum sound power reduction. In the optimi- 
zation process, hundreds of calculations are often 
required to determine the desired configurations. 
Therefore, any significant reduction in the number 
of grid points or elements in a numerical analysis 
will reduce the cost of obtaining the desired 
optimization. 
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Beginning in 1974, Baumeister^ - ® developed a 
wave envelope concept which was used to signifi- 
cantly reduce (two orders of magnitude) the number 
of grid points associated with the "steady state" 
solution of high frequency sound propagation in 
ducts. This concept involved a transformation of 
the wave equation into a form whose solution does 
not oscillate in the axial direction. The use of 
the wave envelope theory' drastically cut the com- 
puter costs associated with the optimization of 
multi-segmented liners. Nayfeh and Kaiser®’^ ex _ 
tended the method for sound propagation in non- 
uniform ducts and with sheared flow. Astley and 
Eversman 111 applied the wave envelope approach 
very successfully in finite element duct sound 
transmission studies. Finally, these same authors 11 
made a very significant extension of the wave enve- 
lope concept to describe simultaneously the induct 
propagation of sound in a turbofan nacelle and its 
subsequent far field radiation pattern. Conse- 
quently, numerical techniques can now be employed 
to study sound propagation in both the internal 
and far field regions of a duct provided the sound 
frequency is reasonably low. 

'v In order to further reduce computer storage 
and run times, as an alternate to the "steady- 
state" theories, time dependent numerical solutions 
were developed by Baumeister 12-1 ^ for harmonic 
sound propagation in ducts. The transient formula- 
tion generally uses a time marching solution to the 
wave equation; consequently, the matrix storage re- 
quirements inherent in the "steady state" formula- 
tion are completely eliminated in the time depend- 
ent analysis. Only the solution vectors for 
pressures and velocities need be stored. By elimi- 
nating the large matrix storage requirements, 
numerical calculations for higher frequency sound 
are now possible. 

The time dependent theory has also been 
applied to forms of the inhomogeneous wave equa- 
tion by Maestrello, Bayliss and Turkel 15 and 
Baumei ster 1 ”. White 1 ' has extended the transient 
theory by means of a mapping to variable area 
ducts. The simultaneous calculation of induct and 
the far field associated with a turbofan engine 
and an unflanged cylindrical duct has been calcu- 
lated by White 1 * 1 and Hariharan and Bayliss 1 ’, 
respectively. 

Considering the grid point reduction of the 
previously discussed wave envelope theory and the 
elimination of matrix storage by the transient 
solution formulation, a logical extension would be 
to combine both theories. In principal, the com- 
puter storage requirements could now be reduced 
many orders of magnitude over previous theories 
making possible calculations with higher frequen- 
cies of three-dimensional fields. 

To combine the theories, the wave envelope 
transient wave equation and appropriate boundary 
conditions will first be derived. The theories 
will be presented for a two-dimensional soft-wall 
duct without mean flow. Next, the complete set of 
difference equations and stability requirements 
will.be presented. Finally, sample calculations 
are presented for plane wave propagation in a hard- 
walled straight duct and for the attenuation pro- 
vided by a two-dimensional soft-walled duct. 


Grid Point Problem 


The propagation of sound in a duct is 
described by the wave equation and appropriate 
source and impedance boundary conditions. The 
wave equation in a two-dimensional rectangular 
duct can be expressed in dimensionless form as: 


1 

C 


2 




2 

ay' 


( 1 ) 


The "steady state" version of this wave equation 
assumes P is proportional to e 1 “ t such that the 
left-hand side is replaced by (iw/c) 2 p. There- 
fore, as mentioned in the introduction, the "steady 
state" form of Eq. (1) is independent of time. 

The usual notation for pressure, distance 
coordinates and speed of sound are used. These 
and other symbols are defined in the nomenclature. 
Here, the dimensionless speed of sound C and the 
dimensionless frequency n are defined as 


1 C* 

_ n " H*f* 


( 2 ) 


The asterisks denote dimensional quantities. 


In a finite difference numerical analysis of 
the wave equation, the continuous acoustic field 
is lumped into a series of grid points (Fig. 1). 
Next, the wave equation is expressed in difference 
form. The difference equations are then solved by 
a time marching process to obtain the pressure at 
each grid point. Obviously, the more grid points 
required in a solution the greater the computer 
storage and run times. 

Consider a hard-wall duct, infinite in extent 
with a plane pressure wave sinusoidal in time at 
x = 0 (P(o,t) = e 11 "* 1 * = e l2l,1 l. The solution of 
Eq. 1 for the pressure yields 12 , Eq. 27 is 

P(x,t) „-i2irr|X ,, 

TZ^T =e . 

e 


The finite-difference approximation _f or the "steady 
state" acoustic pressure in a semi-infinite hard- 
wall duct 12 and the analytical solution (Eq. (3)) 
are presented in Fig. 2 for x greater than 0 and 
less than 1. By a series of numerical experiments, 
the number of axial grid points, I, necessary to 
obtain pressure profiles, velocity and intensities 
accurate to about 4 percent were determined to be 

I = 12 n (L*/H*) + 1 (4) 

Thus, for the case of unit frequency (n = 1) and 
unit duct length to height ratio (L*/H* = 1) shown 
in Fig. 2, 13 grid points were necessary to de- 
scribe adequately the sinusoidal form of the spa- 
tial pressure dependence. If the frequency or 
length is doubled, clearly twice as many points 
will be required to describe the wave, since two 
wave lengths of sound must be resolved. The total 
grid points will be the product of I and the num- 
ber of grid points J in the transverse direction. 

Next, a concept is considered which eliminates 
the direct dependency of the number of axial grid 
points on n and L*/H*. 
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Wave Envelope Equation 


Consider the example case of a soft-wall duct 
with an L*/H* of 3 and an inlet plane wave with 
a dimensionless frequency n equal to 1. A typi- 
cal pressure (real part) profile in a suppressor 
duct is shown in Fig. 3 by the heavy solid line, 
while the dashed lines represent the envelope of 
the solid pressure wave amplitude. The pressure 
amplitude decreases in the axial direction down 
the length of the duct because of acoustic energy 
dissipation at the suppressor wall. However, the 
basic axial variation of the pressure is similar 
to that in a hard-wall duct. Thus, the number of 
grid points (open circles in Fig. 3) still depends 
on frequency n and the dimensionless duct length 
L*/H* (Eq. (4)). 

From Eq. (4), the number of grid points I 
needed would be about 37. However, if the basic 
wave equation could be transformed so that it would 
describe the envelope (dashed lines) of the pres- 
sure rather than the pressure wave itself, the 37 
point requirement could be greatly reduced. As 
shown in Fig. 3, five points are used to describe 
the pressure envelope. 

The assumption is now made that the pressure 
P(x,y,t) can be separated into the following 

P(x,y,t) = p(x,y,t) e“ l2l,n X (5) 

where p(x,y,t) represents the envelope of the 
pressure wave as shown by the dashed line in Fig. 3 
and where 

n + = H*/x* (6) 


some higher order mode, n + would be assumed to be 
a value associated with that mode. 

Governing Equations and Boundary Conditions 

Besides the wave envelope Eq. (7), the equa- 
tions for the acoustic velocity, soft-wall boundary 
conditions, and acoustic intensity are required to 
obtain an expression for the attenuation of a soft- 
wall duct. These equations are now presented. 

Linearized Momentum Equation 

In the absence of mean flow, the x and y 
dimensionless momentum equations can be written as 


3 U _ 
3t 

1 3P 
n 3 x 

(8) 

aV _ 

1 3P 

(9) 

at 

' n ay 


or in terms of the wave envelope parameters 


3u 1 [ap o + 1 

It = - n [ax “ l2,ni P J 

(10) 

3 V 1 3p 

at n ay 

(ID 

Wall Boundary Condition 


In the transverse direction, the acoustic 
impedance at the wall shown in Fig. 1 is defined 
as the ratio of pressure to the transverse velocity 

z* P 

c = p *c* = V 

(12) 


with x* representing the effective axial wave- 
length of the pressure in the duct. Substituting 
Eq. (5) into the wave Eq. (1) yields a new time 
dependent governing differential equation called 
the time dependent wave envelope equation: 


Substituting Eq. (12) into Eq. (9) to eliminate V 
yields 


_aP _ _ n jsP 
ay _ t at 


(13) 


at 




ax 


ay 


+ 

i4nn ap 
ax 



or in terms of the wave envelope parameters 

id — H id (141 

ay " c at ( 1 


In free filed applications, the axial wave 
length x* is known. Thus, the use of the wave 
envelope concept is relatively straightforward. 
However, in soft-wall ducts where multi-modal prop- 
agation occurs, the axial wave length x* is not 
known precisely; therefore, the problem of picking 
n to exactly define the wave length must be con- 
sidered. In the "steady-state" wave envelope solu- 
tion presented 5 , x* in the soft-wall duct was 
assumed equal to x* for a plane wave in a hard- 
wall duct. This assumption was found to yield 
excellent results in the "steady-state" wave length 
solution. Therefore, the same assumption will be 
used here in the solution of Eq. (7). 

As shown in Ref. 5 (Fig. 4), it is necessary 
only to pick a value of x* or (n + ) in the vi- , 
cinity of the average wave length to get large 
savings in the grid points required for a finite- 
difference analysis. Therefore, for a plane wave . 
source as considered herein, n + will be assumed to 
equal n. In a problem where the source might be 


Equation (14) sets the pressure gradient along the 
upper and lower walls. 


Exit Impedance 


In a manner similar to the wall impedance, the 
axial impedance at the. duct exit can be defined as 


e 


P (L*/H*, y,t) 
e - U ( L*7TT* , y,t) 


(15) 


Again, substituting Eq. (15) into Eq. (8) yields 


_3P n aP 

ax ~ c e 3 1 


(16) 


or substituting Eq. (5) into Eq. (16) yields 


3£ = n_ i£ + 

3X q e 3 1 


- A + 

i2nn P 


(17) 
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For the plane wave propagation to be considered 
herein, c e is taken as 1, which is exact for plane 



wave propagation in an infinite hard-wall duct. 

Also choosing to 1 has lead to close agree- 
ment between numerical and analytical results for 
plane wave propagation into a soft-wall duct^. 

Entrance Condition 


The boundary condition at the source plane 
P(o,y,t) allows for transverse variations in pres- 
sure; however, as mentioned earlier, the numerical 
technique will be compared later to previous solu- 
tions in which the pressure and acoustic velocities 
were assumed to be plane waves at the entrance and 
to vary as e 1 ^. Therefore, the source boundary 
condition used here is 


P(o,y,t) = p(o,y,t) = e l2,rt t > o (18) 
Initial Condition 


represented by 
time and space 



usual central differences in 




- (2nn /n) 2 p k j 


( 22 ) 


For times equal or less than zero, the duct is 
assumed quiescent, that is, the acoustic pressure 
and velocities are taken to be zero. For times 
greater than zero, the application of the noise 
source (e.g., Eq. (18)) will drive the pressures 
in the duct. 


where i and j denote the space indices, k the 
time index, and ax, ay and at are the space and 
time mesh spacing, respectively. All spacings are 
assumed constant. The time is defined as 

t k+1 = t k + at (23) 


Acoustic Intensity 

The sound power that leaves a duct and reaches 
the far field is related to the time-averaged 
acoustic intensity, defined as 

I x = 1/2RE |(P/e l2lrt ) x (U/e l2l,t )j (19) 

where the bar represents the complex conjugate. 

The total dimensionless acoustic power is the 
integral of the intensity across the test section 



Solving Eq. (24) for P^ . yields 

■ 9 J 



By definition, the sound attenuation can be written 
as 

adB = 10 log 10 (Ex/Eo) (21) 

Next, the difference form of these equations will 
be presented. 

Difference Equations 

Instead of a continuous solution for pressure 
in space and time, the finite-difference approxima- 
tions will determine the pressure at isolated grid 
points in space as shown in Fig. 1 and at discrete 
time steps At. Starting from the known initial 
conditions at t = 0 and the boundary conditions, 
the finite difference algorithm will march out the 
solution to later times. The development for the 
algorithm that applies to each cell in Fig. 1 will 
now be discussed. 

Central Region (Cell 1) 

Away from the duct boundaries, in cell 1 of 
Fig. 1, the first and second derivatives in the 
wave envelope equation (e.g., Eq. (7)) can be 


a = At 2 / n 2 ay 2 (25) 

6 = 4un At 2 / n 2 (26) 

+ 2 

y = (2uAtn In) (27) 

Equation (24) is an algorithm which permits 
marching-out solutions from known values of pres- 
sure at times associated with k and k-1. The 
procedure is explicit since all the past values 
of P* are known as the new values of k + 1 are 
computed. For the special case at t = o, the 
values of the pressure associated with the k-1 
value are zero from the assumed initial condition. 

Boundary Condition (Cells 2 to 6) 

The expression for the difference equations at 
the wall boundaries are complicated by the impe- 
dance conditions and the change in geometry of the 
cells 2 to 6 in Fig. 1. The governing difference 
equations can be developed by an integration proc- 
ess in which the wave envelope equation (e.g., 

Eq. (9)) is integrated over the area of the cells 
and time. The procedure for the temporal and 
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spatial integration over the cell area is fully 

documented^* 12 . 

. For ease in bookkeeping the solution for 
P K . for the various cells is written in the 
gjfiiral form 


pk+1 a 

P i,j " H - «V 


i 6 k, 




ml k 

p i-l. j 


+ b P k . . + (c + - - i + ij ) p k . 

m 1 ,J-1 I IJ a a a ml ,J 


+ . k + ( 2!jn\ n k 

d m p i,j+l rm o ) p i+l,j 


t 1 -^] k _! 

n^r p ij 


(28) 


condition and are now functions of e^ 2lTt . Thus, 
Eq. (10) becomes • 


i sp , n 
U = 2^ + “ P 


or 


/ k „k 

(Pi+i h ~ Pi_ 


u‘ , = ± i - 1>J ) * 4- p" 


__ lzLi]+ ,+ 

J i ,j ” 2nn y 2ax J n r i,j 

The intensity, (Eq. (19)), can be written as 


T k 1 D -k . i 2irt k , i2irt 

I-J^RePi/ xu i,j /e 


where the bar over p represents the complex 
conjugate. 


(31) 


(32) 


(33) 


where the values of the coefficients a™, b m , etc. 
for the various cells are listed in Table I. In 
particular, substituting the values of the coeffi- 
cients for cell 1 in Table I into Eq. (28) will 
reproduce Eq. (24). 

Stabi 1 ity 

In the explicit time marching approach used 
here, round-off errors can grow in an unbounded 
fashion and destroy the solution if the time incre- 
ment at is taken too large. The von Neumann 
method is often used to study the stability of the 
difference approximations to the wave equation. 
Application of the von Neumann method 21 (pg. 104) 
to Eq. (22) requires that the time increment be of 
the form 



for n + small, Eq. (29) reduces to the conven- 
tional CFL (Courant, Friedrichs and Levy) condi- 
tion. Because of possible effects of boundary con- 
ditions, some numerical experimentation was used to 
check the validity of Eq. (29). 


Next, the total dimensionless power is deter- 
mined according to Eq. (20) 

/ LAST-1 \ 

E ( ■ (t 'll * Z ‘5.J * « <»> 


and finally the attenuation is 
adB = 10 log 


(E k> 


lOlrk 
1; 


(35) 


Sample Calculations 


In two sample problems to follow, the time- 
dependent wave envelope results will be compared 
to closed form analytical' solutions. First, the 
simple case of plane waves propagating down a hard- 
wall one-dimensional duct is presented. This case 
allows comparison of the numerical and analytical 
pressure profiles down the length of the duct. The 
second example compares the numerical and analyti- 
cal predictions of the attenuation in a soft-wall 
two-dimensional duct. 

One-Dimensional Hard-Wall Duct 


"Steady State" 

Recall, at the start of the numerical calcula- 
tion, the acoustic pressures and velocities were 
assumed zero throughout the duct and a pressure 
source begins a harmonic oscillation at x = o for 
t > o. The transient numerical calculation must 
be started and continued in time until the initial 
acoustic transient has died out. This subject has 
been treated extensi vely 1 ^’ . For plane wave 
propagation, as shown in these references, the 
transient solution for pressure equals the "steady- 
state" Fourier transform solution when 

t > n (L*/H*) (30) 


Numerical and analytical values of the pres- 
sure were computed for the special case of a hard- 
wall duct with a L*/H* of 1, n of 1 and an 
inlet plane wave. The analytical value from 
Eq. (5) yields 

= 1 + io *o 06) 

e 

As seen in Fig. 4, the agreement between the 
numerical and analytical results is reasonable. A 
comparison with Fig. 2 indicates the essential dif- 
ferences between the transformed numerical solution 
and the conventional numerical solution for the 
same problem. 


Axial Velocity and Intensity 

The axial velocity (Eq. (10)), acoustic in- 
tensity (Eq. (19)), and the sound attenuation 
(Eq. (21)1 can be quite easily expressed in differ- 
ence fornr - °. In these cases, the pressure and 
velocity are now assumed to be at "steady state" 


Two-Dimensional Soft-Wall Duct 

As another example of the transient wave 
envelope formulation, the noise attenuation at the 
optimum point (point of maximum attenuation in the 
impedance plane) is now calculated for a two- 
dimensional duct with L*/H* values varying 
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between 0.5 and 6 and input plane waves with dimen- 
sionless frequencies n of 1, 2, and 5. This 
range of dimensionless parameters essentially 
covers the practical range of application to turbo- 
jet exhaust suppressors. 

For a plane entrance pressure profile the 
closed form analytical results show that the values 
of specific acoustic wall impedance c associated 
with the optimum impedance used in the numerical 
analysis were determined from the analytical 
techniques 23 ’ 24 and are listed in Table II. 

The numerically calculated attenuations are 
compared to the corresponding analytical 
results 23 * 24 which are applicable to infinite 
ducts. The numerical results (symbols) and the 
analytical results (lines) for the maximum attenu- 
ation are shown in Fig. 5. The analytical and 
numerical results are in very good agreement. A 
slight deviation occurred at the very high attenu- 
ation associated with the low frequency case n = 1 
at L*/H* of 4. For the low frequency case, it is 
suggested that more points be employed. 

Based on the wave envelope concept, the numer- 
ical calculations shown in Fig. 5 used 30 grid 
points in the axial direction. For n = 5 and 
L*/H* = 6, the standard finite difference technique 
required 360 grid points, according to Eq. (5). 
Thus, for a J = 10, the total number of grid 
points has been reduced from 3600 to 300 when the 
wave envelope concept is employed. This represents 
an order of magnitude savings in computer storage 
and computational time. 

Concluding Remarks 

Transient finite difference solutions using 
the wave envelope concept are presented for plane 
wave sound propagation in a one-dimensional hard- 
wall duct and a two-dimensional soft-wall duct for 
zero Mach number. The results show the numerical 
procedure to be in agreement with the corresponding 
exact analytical results. 

The wave envelope approach to the numerical 
problem reduces the number of grid points in the 
difference solution by an order of magnitude com- 
pared to the conventional difference technique. 
Table III shows the large reduction in computer 
storage requirements by employing both the tran- 
sient and wave envelope techniques as compared to 
the conventional "steady-state" theory. Clearly, 
numerical solutions for acoustic propagation in 
complex ducts can now be determined with reasonable 
storage and computational times. 
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TABLE II. - OPTIMUM IMPEDANCE 
AND ATTENUATION FOR UNIFORM 
RECTANGULAR INFINITE DUCT 
















































TABLE III. - GRID POINT AND STORAGE REQUIREMENTS (REAL 
AND IMAGINARY) FOR n = 6, L*/H* = 5, J = 10 
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Figure 1. - Grid-point representation of two dimentional duct 



Figure 2. - Analytical and numerical pressure 
profiles for one dimensional sound propa- 
gation in hard-wall duct, n= 1. 
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Figure 3. - Typical pressure profile for sound prop- 
agation in a soft-wall duct for dimensionless 
frequency n = 1 and duct length L*/Ff = 3. 
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Figure 4. - Analytical and numerical dimension- 
less acoustic pressure profiles for sound 
propagating in a hard wall duct. (r| = 1, 

L*/FT = 1, 1 = 5, J - 10, t= 5.0, At = 

0. 5 x At max ). 
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Figure 5. - Effect of axial length and 
frequency on attenuation at optimum 
impedance in two-dimensional duct 
for a plane wave input (I = 31, 

J ■ 20, t-5.0, At = 0.5xAt max ). 
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